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Abstract 

We report analytical calculations for the propagation of electromagnetic radiation through an inhomogeneous layer 
whose refractive index varies in one dimension situated between bulk right- and left-handed media. Significant field 
localization is generated in the layer that is caused by the coherent superposition of evanescent waves. The strength 
of the field localization and the transmission properties of the layer are investigated as a function of the layer width, 
losses and defects in the refractive index; the former two being modelled by continuous changes, and the latter by 
discontinuous changes, in the index profile. 
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1. Introduction 

Metamaterials are of importance because of their ability to manipulate the behaviour of electromagnetic 
radiation in ways that cannot be achieved by naturally occurring matter |1I2I3| . They function through 
fabricating the microscopic electromagnetic properties to create an effective medium, and Maxwell's macro- 
scopic equations provide the description of the fields therein. A class of metamaterials with permittivity and 
permeability simultaneously less than zero have been fabricated [41516] , these having been postulated earlier 
[7]. These left handed media (LHM) are so called because their E, H and k vectors form a left-handed 
co-ordinate set with the consequence that their wave propagation vector is anti-parallel to the direction of 
energy flow. The unusual optical properties of such materials have stimulated great interest, not least be- 
cause of their ability to create lenses that are not diffraction limited, and for cloaking T"2'8], whilst providing 
a striking example of resonant behaviour in physics However, the resonant behaviour implies that the 
interaction of radiation with LHM is a singular phenomenon and is therefore susceptible to imperfections in 
the material parameters and the way that they are configured in bulk matter. 

Left handed media are necessarily dispersive [9] which leads to losses that affect substantially the efficiency 
of components constructed from them. The values of the material parameters required for their applications 
to operate successfully are exacting |10lll| . and deformations and contamination of the idealized planar 
surfaces between components can also severely limit performance |12I13| . This paper examines the effect 
of an imperfect boundary between bulk right- and left-handed media, such as may be caused by a slight 
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Fig. 1. (Colour Online) (a) A typical one-dimensional GRIN layer profile containing a step discontinuity of 5. (b) Pictorial 
representation of the incoming planar wave and the resulting reflected and transmitted waves. 

roughening of the interface, or due to homogenization of the material parameters through the interminghng 
of bulk material of different effective refractive index. 

The model adopted will assume that the thickness of the transition layer between the bulk right- and 
left-handed material is small compared with the smallest lateral scale-size that characterizes the interface. 
Consequently the change in refractive index can be assumed to be limited to one-dimension. Maxwell's 
equations then reduce to an ordinary differential equation, for which the solution is expressible in terms 
of exponential functions in the bulk media, and in terms of special functions for particular profiles of the 
refractive index within the layer. In particular a linear transition between refractive index of -1-1 and -1 
has an analytical solution within the layer in terms of confluent hypergeometric functions, and it is the 
properties of this exact solution that are explored in this paper. 

The following section formulates the model for the layer and discusses the general forms of the solution 
to Maxwell's equations. Section [3] obtains the solution for the linear refractive index profile and discusses 
both the form and interpretation of the solution in the absence of losses in the material parameters. The 
most striking feature of this solution is a localization of the field within the layer. Section |4] reconsiders the 
model in the presence of finite losses and shows that the lossless solution is singular in several respects, 
not least that losses destroy the localization and perfect transmission property of the layer. Section O shows 
that a small discontinuity in the refractive index profile obtains solutions whose properties are intermediate 
between the lossless and lossy cases. The concluding section summarizes and discusses the results. 

2. The GRIN model 

Graded-index (GRIN) modelling of a boundary between bulk media has been used in conventional RHM 
for many years 14J and is tantamount to considering the properties of waves that interact with a one- 
dimensional refractive index profile 



where i^izi) = ni, iy'{z2) — «2 with m and n2 being the bulk refractive indices in media 1 and 2 respectively ^ I. 
The refractive index in the layer can be continuous, piecewise continuous or contain discontinuities of size 6, 
as shown in figure [1] (a). There is little restriction on the profile of viz), but only a few analytical solutions 
of Maxwell's equations with profile ^ exist for problems involving LHM. 

^ More speciflcally e{zi) = ei, e(z2) = £2, /^(zi) = Ml and /^(z2) = with subscript 1 indicating medium 1 and subscript 2, 
medium 2 




n2 Z2 < Z 



v{z) Zi < Z < Z2 



Z < Zl 



(1) 




2 



A harmonic s-polarized electromagnetic wave travels in the direction of increasing z through a lossless 
right handed medium of refractive index ni, and at angle B to the normal of the transitional layer, as shown 
in figure [1] (b). The form of the electric field is E(a;,y, z) = Ey{x,y, z)y (a harmonic time dependence, u> 
is assumed throughout). Because the refractive index is a function of z only. Maxwell's equations lead to 
separable solutions of the form Ey(x, y, z) = exjp{ikxx) E(z), where is the wave number in the x direction 
which is in the plane of the interface between the bulk media. In the homogeneous bulk media the form 
adopted by E{z) are simple exponential functions, whereas in the layer the field is a solution of the equation 

which results in the electric field adopting the form 

exp(i /c^j z) + r exp(— z k^-^ z) z < Zi 
E{z) = { E{z) z^<z<Z2 , (3) 

t exp(i kz2 z) Z2 < z 

where E is the electric field in the layer, c is the speed of light in vacuum, fc^^ and k^^ are the wave number 
in the z direction in media 1 and 2 respectively. The reflection and transmission coefficients r and t are 
determined once the form for E is known by invoking the continuity of the tangential components of E and 
H at both extremities of the layer at zi and Z2. 

A number of profiles for e{z) and ^i{z) give rise to analytic solutions to equation Of these, there are 
two that can traverse between RHM and LHM: that of an exponential profile and that of a linear profile. 
Only the linear profile allows for independent choice of the losses within the model, and therefore will be 
adopted. The results of the exponential case are given in Appendix [Al The profiles for ^{z) and e(z) for the 
case of interest are taken to be 

fi(z) = m, z + d, e{z) = T] fi(z) . (4) 

It is clear from these definitions that with complex m and d, this profile can be used to model the transition 
between two bulk media of any /i, whether lossy or not. The limitation that e cannot be chosen independently 
of /i is not especially restrictive given that most devices which employ negative refraction or perfect lensing 
have e and /i everywhere equal [118] . i.e. that rj is unity in equation ([¥]). Consequently rj will be set to unity 
in the remaining sections. 

For the profile given by ^ the electric field within the layer is given by 



exp(^i^^)vl'(^) 

where 



^(^) = ,22 ("^(^) + /^GC^)) (5) 



F(2) = M( 1-^^,2,17^(2) 
4m"^ 7 

G{2) = U( l--^,2,i7*(2) 
\ 4m'' 7 

„i/2 ^ 

7 = , <!f{z) = (d + mzf = fi'^, 

cm 

M and U are the independent confluent hypergeometric functions [T5^ and a and (3 are constants of 
integration. Using the conditions that 

Eiz) and (6) 
/i(z) dz 

are continuous across both boundaries between the inhomogeneous layer and the bulk media uniquely de- 
termines the four unknowns, r, i, a and (3. These have a similar structure to those given in [14J. In the 
first instance this will be used to determine the electric field throughout a diffuse boundary between lossless 
RHM and lossless LHM. 
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Fig. 2. (Colour Online) The real part of E that results for a dimensionless layer width, ka = 8tt for an angle of incidence 
9 = 7r/8 to the normal of the surface from a vacuum (z < —1) to a bulk medium with e = /i = —1 (z > 1). (a) shows a 3D 
representation of the real part of the electric field, note the localizations move with speed uj/kx in the x direction, (b) the same 
situation in contour form and (c) a cross section of part (a) in the constant x plane. The amplitude of the field is taken to be 
unity in the homogeneous media. 

3. The lossless solution 



Without loss of generality, the layer is assumed to be of width 2a and is centred at the origin. Figure [2] 
shows various depictions of the real part of the electric field when the s-polarized wave is incident from a 
right-handed medium (n = 1) into a left-handed media {n = —1) at an angle of incidence 6 = tt/S, and 
when the dimensionless layer width ka = 8tt. 

The surface plot figure [2] (a) illustrates the localized field confined within the transition layer. Figure [2] 
(b) depicts the level contours of the field and clearly shows that negative refraction has occurred from the 
orientation of the phase fronts far from the layer. Figure [2] (c) is a cross-section across the layer, and shows 
that the field within the layer is enhanced several times in excess of its (arbitrary) value in the homogeneous 
media. The localized structure is, moreover, symmetrical about the origin, there being no way that an 
asymmetrical solution can match onto the field structures in the bulk media. 

Although the refractive index changes smoothly from +1 to —1 with increasing z, the impedance Z = 1 
everywhere, and so the layer is perfectly transparent to the radiation, i.e. t = 1 and r = 0. Consequently if 
the fields were observed in the bulk media alone, it would not be possible to infer the presence or absence 
of the layer. 

The origin of the localized field can be understood with reference to figure 3, which is meant for illustrative 
purposes alone and which depicts the field through the layer, within which \n\ < 1. In the region — l<z/a< 
0, both e and fi are positive and so when the value of |n| falls below that which allows a propagating mode 



4 




Fig. 3. (Colour Online) Within the inhomogeneity layer there is a band where an incoming wave becomes evanescent in nature, 
\n\ being too low to allow propagation (a). These evanescent modes decay within the RHM and exponentially grow in the LHM 
[8] (b). When an evanescent wave reaches a point where it can propagate again the boundary conditions cause an evanescent 
wave to be reflected back into the layer (c). This propagates to the other side where another reflection occurs (d). It is the 
summation of each of these evanescent reflections that gives rise to the structure seen in flgure|2](c). 

to exist, i.e. \n(z)\ < sin0, an evanescent mode is established that decays with increasing z. This evanescent 
mode has a finite amphtude at z = 0. For < z/a < 1, both e and /i are negative, with result that 
the evanescent wave is amplified out to a distance Zc satisfied by |n(2:c)| = sinfl, whereupon the wave can 
propagate once again for z > Zc- An evanescent wave is refiected back into the layer that increases with 
decreasing z, until z < 0, where the mode then decreases in the right-handed medium. Thus a coherent 
structure is established in the region |?i(z)l < Zc contained within the layer through the interference of these 
evanescent modes. This prompts questioning how the peak magnitude of the field scales with the layer width. 
This dependence is shown in figure 2] (a) by the red (triangles) curve, and shows an exponential growth with 
the layer thickness. This dependence obtains from the complicated nature of the confluent hypergeometric 
functions and their derivatives, that are contained within a and /3. Clearly this behaviour is unphysical 
and will be ameliorated by the incorporation of losses, as will be considered in the next section. The red 
(triangles) curve in figure [IJb) shows the peak value of the field as a function of the angle of incidence for a 
dimensionless layer width ka — Svr. This increases with increasing 9 since the optical distance through the 
amplifying layer increases with increasing layer width. 

4. Lossy Transition Prom n — 1 to n — —1 + Ki 



In this section losses are incorporated into the bulk left-handed medium, which is assumed to have refrac 
five index n — —1 -I- ni with k > 0. The inhomogeneous layer is also assumed to be lossy with profile: 



fi{z)^e{z) = + £ + 



(7) 



According to this model, the refractive index is purely imaginary at the origin and the material is no longer 
impedance matched for all values of z. These two properties imply that the field will no longer grow without 
bound with increasing layer width and that the layer is no longer perfectly transparent to the radiation, 
i.e. a reflected wave exists in the right-handed half-space, z < 0. Indeed the presence of losses implies that 
for z ^ 0, the wave will have decayed completely. The structure of the solution within the layer is also 
substantially modified by the losses, for the second confluent hypergeometric function possesses a branch 
cut which must be crossed as the layer is traversed. It can be shown that for a change in refractive index 
from ni = (f>i + Kii at zi to n2 = (f>2 + «^2* at Z2 then the position z where the branch cut is crossed is given 
by: 

h - K2 </>!) |Az| + (Ak2 + A(/.2)'/' (zi (/)2 - Z2 0l) 



(Ak2 + A02)1/2A(/. 

with A-0 = ip2 — 4'i for any ■0- For z > z, U can be calculated through the use of 



U(a, 2, Z(z)) = U(a, 2, Z{z)) - H(z - z) 



2TTi Sign [Me (n|^^-)] 



M(a, 2, Z(z)) 



(8) 



(9) 



r(a-i) 

where z is given in a and Z{z) are the first and third argument of the U found in equation ([5]), H is 
the Heaviside step function and U is the principal branch of the hypergeometric function [13] that is defined 
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Fig. 4. (Colour Online) A wave of wave number k, propagating at an angle d to the normal, transitions from a vacuum through 
the inhomogeneity region, of width 2a, into a medium where = —1 + ni. (a) shows the magnitude of the localisation in the 
medium relative to the incident wave with 9 = tt/8 for k = (lossless) [Red Triangles] and k = 10~^ [Blue Squares], (b) also 
shows the magnitude of the localisation but as a function of the incident angle, 9, for no losses [Red Triangles] and k = 10~^ 
[Blue Squares] 
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Fig. 5. (Colour Online) For the situation detailed in figure [J] (a) displays the dependence of the transmission coefficient on 
9 for various ka, this is obtained with k = 10~^ for ka = 5 [Red Squares], 1 [Blue Circles] and 0.1 [Black Triangles] and (b) 
shows the transmission coefficient as a function of fca at 6 = n/S for k = 10~^ [Red Squares], 10""^ [Blue Circles] and 10"^'^ 
[Black Triangles] (The latter two are indistinguishable) 

when — TT < argZ(z) < tt. This allows for calculations involving losses to be performed in a similar way to 
the previous sections. Figure [4] details a variety of different features resulting from the addition of losses. 
The main point to note is that figure [¥] (a) shows the exponential dependence of the localisations on layer 
width is suppressed by the addition of losses. 

A careful treatment of the branch-cut shows that the form of the solution is affected. The blue (squares) 
curve in figure 3] (a) shows that the peak value of the field in the layer now decreases exponentially with 
increasing layer width for a modest value of k = 10^"^, as is also the case for the variation with 9 as shown 
in figured] (b). Indeed the peak value of field within the layer is less than that outside it, so that the 
localization effect is entirely suppressed. This is because the evanescent modes are dissipated by the losses. 
Despite these losses, the field within the layer is still finite once the location where the modes can propagate 
again is attained, and so there is still transmission of radiation into the bulk left-handed medium. Figure 
[5] (a) shows the transmission coefficient at z = a for k = 10^^ as a function of the angle of incidence for a 
selection of values of the dimensionless layer width, ka = 0.1 (brown triangles), ka — I (blue circles) and 
ka — 5 (red squares). The layer becomes opaque to radiation at progressively smaller angles of incidence as 
the layer thickness increases. Figure [5] (b) shows the negative exponential dependence of the transmission 
coefficient as a function of ka for different values of the loss, k = 10^^ (red squares), k = 10^'' (blue circles) 
and K = 10^^" (blue triangles), the latter two being essentially indistinguishable. Hence for losses of the 
order 10^'^ or lower the transmission coefficient is insensitive to the precise value of k: this is an artefact 
of the dampening of the evanescent reflections within the transition layer. It can also be noted that the 
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transmission coefficient, for fixed ka, increases as the losses increase. Tfiis effect occurs because large losses 
quickly damp all but the first reflection of the evanescent mode within the layer which leads to a larger ratio 
between the transmitted and reflected wave. The total transmitted and reflected power, however, steadily 
decrease as losses increase. 



5. The effect of discontinuities in the refractive index profile - the ^-GRIN model 

This previous section showed that the lossless results are very different from those where losses are included 
and this is principally because the refractive index vanishes at z = in the former case, whereas there is 
a branch-cut in the solution for the latter. This section will examine the robustness of these two classes of 
solution by modelling the refractive index in the layer by a staircase. Within each plateau of the staircase, 
the refractive index is constant and so the solution to Maxwell's equations is comprised of two independent 
exponentials. Using the boundary conditions ([6|) at the end of each of the N steps gives 2N equations for 
the 2N constants that determine the amplitudes and phases throughout the layer [l4] . 

Figure [6] contrasts the solutions obtained from the continuum models of sections [2l|4] and the discrete 
staircase model. Figure [S^a) is for a lossless medium and is repeated from figure [^c) for ease of comparison; 
this should be compared with[6Ub), which is for a staircase with N = 2^ equally spaced plateaus throughout 
the layer, and there is no discernible difference between the two solutions. Figure ^ic) shows the GRIN 
model solution for k = 10~^ - the localization has been suppressed entirely by the losses. Note however that 
the solution of the lossy staircase model shown in[n)Jd) retains the localization feature and is essentially an 
attenuated form of the solution shown inlHIb). Hence the staircase model is quantifiably different from the 
lossy version of the GRIN model. This prompts investigating whether a simple element can be incorporated 
into the GRIN model that captures both the lossless and lossy behaviours shown by the staircase model. 

The element that makes both GRIN and staircase models qualitatively consistent with each other is the 
inclusion of a (single) discontinuity of size 5 located at z = z, see equation ([8|), as depicted in figure [TJa). 
Figure [7] quantifies the sensitivity of this (5-GRIN model with the size of 5 for an angle of incidence — 7r/8. 
Figure [7fa) shows the dependence of the size of the peak value of the localized field with the dimensionless 
layer width. The blue (stars) and green (diamonds) curves are the lossless and lossy [n = 10~^) cases 
respectively, as previously seen in figure UJa), and these act as bounds for the lossy (5-GRIN model, where 
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Fig. 7. (Colour Online) (a) Localisation height with the parameter values of figure|4](a) but with the addition of a discontinuity, 
5. The case shown are k = 10-3 and (S = [Green Diamonds], k = 10'^ with <5 = 0.01 [Blue Circles], k = 10~^ with 5 = 0.05 
[Purple Triangles], ft = 10" with S = 0.1 [Dark Yellow Squares] and k = with no S [Light Blue Stars], (b) Transmission as 
a function of S for ka = 1 [Purple Triangles] , ka = 10~^ [Dark Yellow Squares] and ka = 10^^ [Blue Circles]. 



tlie value oi S — 0.1 dark yellow (squares), S = 0.05 purple (triangles) and 6 = 0.01 blue (circles). It can 
be seen that these results match smoothly to the lossless results for small values of ka, and for sufficiently 
small values of ka are independent of S. Figure [71(b) shows the dependence of the transmission coefficient as 
a function of S (note, a value of (5 = 2 is equivalent to step change between the bulk right- and left-handed 
media without a diffuse layer). The curves displayed are for different values of the dimensionless layer width, 
ka = 1 purple (triangles), ka = 10~^ dark yellow (squares) and ka = 10~^ green (circles). Thus the more 
diffuse the layer is, the less radiation is transmitted to the left-handed medium. 



6. The simultaneous limit S — *■ 0+ and k — > 0+ 

We have shown that the GRIN model displays different behaviour according to whether k = or is 
finite, but that the latter case can be made consistent with the staircase model with the inclusion of a 
discontinuity S in the refractive index profile, and this we have termed the (5-GRIN model. There is clearly 
singular behaviour in the nature of the solutions obtained with these models as 6 and k both tend to zero, 
which recovers the lossless GRIN model. This section discusses the reason for the discrepancy. 

We may define two classes of solution obtained in terms of whether or not there exists localization of the 
field within the layer, and the class that the solution adopts is different according to whether we take S ^ 
followed by K — > (no localization), or k ^ followed by 5 (localization). If 5 ^ k (in magnitude) as 
(5 — > 0, then the field will show pronounced localization features. 

The terminology of Lconhardt 16J provides an alternative way to view the behaviour of the two models, 
depending on whether the branch cut of the U hypergeometric function is crossed. The quantity that 
determines this is the third argument of U, which is A(z) — i'j'^iz), see equation ([5]). In the lossless case 
the path of A(z) causes the third argument of U to touch its branch cut at the origin, but not necessarily 
cross it, this gives rise to the two distinct solutions seen previously. The axes of figure [8l[a-b) are the real 
and imaginary parts of A(z) and the curves are the loci of A(z). The branch cut of U is Re(A(z)) < and 
Im(A(2:)) = 0, which is shown by the thick black lines. Figure [3^ a) is for the GRIN model, and the paths 
of A(z) must cross the branch-cut, so that the solutions move to the next Riemann sheet, and give rise to 
the non-localized solutions as depicted in figure [5](c). Figure [5Kb) is appropriate for the (5-GRIN model, and 
the presence of a non-zero value of S means that the locus of A(z) does not cross the branch cut, but can 
'jump' across it with the solution remaining on the same Riemann sheet - hence the localization form of the 
solution is obtained as shown in Figure [81(d). 
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Fig. 8. (Colour Online) An illustration of the complex path for the third component of the U as the layer is traversed. The thick 
black line indicates the branch cut and the paths are for different losses, (a) shows the GRIN solution and that as k — > the 
path lines up to the negative imaginary axis whilst including the branch cut contribution, (b) shows that the introduction of 
the discontinuity causes the branch cut to be avoided (dashed line) leading to a lossless solution along the negative imaginary 
axis that does not include the branch cut. The real parts of the electric field across a lossless layer are shown in (c) and (d) for 
the limit of (a) and (b) respectively, using typical values for any variables used 

7. Conclusion and discussion 

In this paper we have introduced a model to account for smooth changes of permittivity and permeability 
across a diffuse boundary between bulk right- and left-handed metamaterials. A full-wave, exact analytical 
solution to this problem leads to a strong localisation of the field in the transition region whilst being fully 
transmissive. In the lossy case the strong localisation is removed and a reflected wave exists. Consideration 
of another analytical model leads to the inclusion of a discontinuity in the refractive index profile which 
restores qualitatively the features of the lossless case. In all cases this paper has been able to analytically 
quantify the reflected and transmitted wave properties. 

The reason for the localization in the layer is the constructive interference of evanescent modes that are 
stimulated whenever |n| < sin6' as illustrated in figure ID It should be stressed that these modes are not 
a conventional plasmon mode which is generated by a discontinuous change in the refractive index [17j . In 
contrast the diffuse layer causes the coherent addition of a plasmon and an anti-plasmon [18j throughout 
the volume of the metamaterial for which \n\ < sm9. 

Although not detailed here, equation ^ can also be used to model changes between two right-handed 
media or indeed two left-handed media. Altogether this method can model a diffuse boundary between any 
combination of left- and right-handed media, with or without losses in either medium. 

Another practical application of the model is the extension of the graded index approach to modelling 
surface roughness as a graded index change, e.g. |14j . to include magnetic and left-handed materials. Also, 
since the solution contains the polarization-state of the wave the approach can be used to investigate, 
for example, the emission polarization effects of infra-red radiation, e.g. jl9j . from left-handed media. Other 
applications include further polarization effects, such as the Brewster angle, and to analyse how superlensing 
properties are sensitive to the polarization state and by roughness of the lens surfaces, and these will be 
treated elsewhere. 

During the period that this paper was being refereed, a paper |20j has studied a similar model to that 
presented here but with different emphasis. 
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Appendix A. Exponential Layer dependance 



The other most useful profile, is that of an exponential dependence across the layer: 

/i = exp(TOz), e — 7]fi + A (A.l) 



with the solution 



where 



Eiz) = (^^) ' exp (-Ziz)) {aF{z) + fiG{z)) (A.2) 



i Ac 



F{z) = V\ V1^^^^2Z(z) 



G{z) = h\ VIJU^^^. 1,2 Z(z) 



1/2 

Ak^\ -r^r \ « exp(m 2)»?^/^aj 



U and L being the Hypergeometric U and Laguerre L functions respectively and a and (3 are constants of 
integration. In the case that A = Q, i.e. that e(z) = r] ii{z) = rj exp(mz) then the U and L functions simplify 
to the Bessel J and Y functions. The most useful applications of this model have /i and e continuous across 
z = zi and z = Z2, c.f. section [21 This allows for a large number of possible changes between two RHM, but 
for a change from a RHM to a LHM m must be complex and necessarily introduces losses into the model 
which are largely uncontrollable. Still this model exists and may be better suited to some applications. 
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